Draft version March 8, 2011 

Preprint typeset using I^TgX style emulateapj v. 08/22/09 



THE DYNAMICS OF SPIRAL ARMS IN PURE STELLAR DISKS 



o 

<N 



in 

< 

6 

> 
oo 

<N 

O 
O 



X 



M. S. Fujii^ J. Baba, T. R. Saitoh, J. Making, and E. Kokubo 

Division of Theoretical Astronomy, National Astronomical Observatory of Japan, 
2-21-1 Osawa, Mitaka, Tokyo, 181-8588 

AND 

K. Wada 

Graduate School of Science and Engineering, Kagoshima University, 
1-21-35, Korimoto, Kagoshima, 890-0065 
Draft version March 8, 2011 

ABSTRACT 

It has been beheved that spiral arms in pure stellar disks, especially the ones spontaneously formed, 
decay in several galactic rotations due to the increase of stellar velocity dispersions. Therefore, some 
cooling mechanism, for example dissipational effects of the interstellar medium, was assumed to be 
necessary to keep the spiral arms. Here we show that stellar disks can maintain spiral features for 
several tens of rotations without the help of cooling, using a series of high-resolution three-dimensional 
A/'-body simulations of pure stellar disks. We found that if the number of particles is sufficiently large, 
e.g., 3 X 10^, multi-arm spirals developed in an isolated disk can survive for more than 10 Gyrs. 
We confirmed that there is a self-regulating mechanism that maintains the amplitude of the spiral 
arms. Spiral arms increase Toomre's Q of the disk, and the heating rate correlates with the squared 
amplitude of the spirals. Since the amplitude itself is limited by Q, this makes the dynamical heating 
less effective in the later phase of evolution. A simple analytical argument suggests that the heating 
is caused by gravitational scattering of stars by spiral arms and that the self-regulating mechanism 
in pure-stellar disks can effectively maintain spiral arms on a cosmological timescale. In the case of a 
smaller number of particles, e.g., 3 x 10^, spiral arms grow faster in the beginning of the simulation 
(while Q is small) and they cause a rapid increase of Q. As a result, the spiral arms become faint in 
several Gyrs. 

Subject headings: galaxies: kinematics and dynamics — galaxies: spiral — methods: n-body simula- 
tions — stellar dynamics 



1. INTRODUCTION 

The physical origin and evolution of spiral arms in 
disk galaxies are long-standing problems of galactic as- 
tronomy. The most widely known theory is the Lin- Shu 
hypothesis, in which spiral structures are interpreted as 
stationary densit y waves with a constant pattern speed 
in a stellar disk (iLin fc Shul ll96^lBertin fc Linl[T996. ). 
However, as pointed out by iToomrd (|1969f ). Lin-Shu's 
mechanism has a serious problem, such that the en- 
ergy and angular momentum of the tightly wound spi- 
ral waves radially propagate with the group velocity, 
and the waves are absorbed at the inner Lindblad res- 
onance. Therefore, a continuous generating mechanism 
is necessary to maintain the stationary density waves, 
e.g., WASER mechanism; an outward-traveling wave is 
reflected and transmitted into other traveling waves at 
the co-rotation resonance (Mar3IT976f ). Alternatively, it 
has been proposed that spiral arms grow from small-scale 
disturbances through t he swing amplification mechanism . 
iGoldreich & Lvnden- Bell 1961; lJulian fc To omrg [19661 : 
pToomrel il981). In this picture, spiral ar ms are recur- 
rent and transient rather than stationary ()Toomrelll99Ql : 
iToomre fc Kalnajl Il991f ). Previous A^-body simula- 
tions supported the recurrent and transient spirals, espe- 



cially for multi - arm spirals (iSellwood fc Car Iberd 119841: 
ISellwood I [200(1 ISjUroo fc BinnevI 120021 : iFuchs et al.l 
l2005l: ^ Sellwoodl l2010f ) and also for barred-spirals 
(|Baba et al.n2009l ). 

While A/'-body simulations both with and with- 
out gas show ed the recurrent and transien t spiral 
arms (Sellwood fc Carlberg 1983; iCarlberg fc Fr eedmanI 
19851: Elmegreen fc Thomasson 1993; BottemgTllQQsF 
Baba et aL .20091 ..Sellwood fc Carlberg. (,.1984 ) pointed 
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out that spiral arms in pure stellar disks (i.e., without 
gas) disappeared in several galactic rotations. They per- 
formed two-dimensional A/'-body simulations and showed 
that stars scattered by spiral arms heated up the disk 
(increase Toomre's Q value), and thereby spiral arms 
disappeared. They argued that some dynamical cooling 
mechanism was necessary to maintain the spiral arms. 
They showed that when new stars with circular orbits 
(i.e., with zero velocity dispersion) were added to the 
stellar disk with a constant rate, the spirals were main- 
tained for about 10 galactic rotations. This demonstra- 
tion was based on the idea that stars are formed from 
the interstellar medium (ISM) with a small velocity dis- 
persion. After their work, the effects of gas and star 
formation wer e investigated for cooling and also heating 
(lElmegreen fc Thomasson 1993; Bottema 2003). In ad- 
dition, 'Bottema (2003) proposed th at the filamen ts of 
gas trigger the swing amplification ()Toomrelll981f ) and 
enhance stellar spiral arms. 
However, the dynamical effect of the ISM to stellar 
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spirals is not well understood yet. Recently, iBaba et all 
(|2009) showed that the multi-phase ISM in a stellar disk, 
in which spiral arms are self-excited, has a complicated 
velocity field, and the cold, dense gas like the giant 
molecular clouds (GMCs) has large non-circular motions. 
Therefore, the newly born stars are not necessarily dy- 
namically cold. Furthermore, since the mass fraction of 
gas in galactic disks is typically only ~10%, it is nat- 
ural to assume that the stella r component controls the 
dynamics of the disk. Indeed, "Elmegreen & Thomasson" 
fl993) performed two-dimensional A^-body simulations 
including gas particles and concluded that the stellar Q- 
value contro ls the formation of sp ir al str uctures. 

Although iSellwood fc Carlberd (|1984f ) reported that 
the stellar disks are heated up significantly in sev- 
eral galactic rotations, we should be careful on the ef- 
fect of numerical artifacts. In particular, the number 
of particles used in their A/'-body simulations is only 
2x10^. In such simulations, two-body relaxation might 
have significantly enhanced the decay of spiral arms. 
iDonner fc ThomassonI (|1994[ ) also performed similar sim- 
ulations for stellar disks with m = 2 spiral arms. The 
number of particles was 5 x 10^. They argued that 
their spiral arms were long-lived. However the lifetime 
was only several galactic rotation s. Their simulations 
have t he same problems as those of lSellwood fc CarlbergI 
(|1984l ). Therefore, three-dimensional A/'-body simula- 
tions with a large number of particles are necessary to 
investigate the long-term evolution of stellar disks. Such 
simulations are now feasible, thanks to the progress of 
computers and numerical methods. However, recent sim- 
ulations of galactic disks ha ve focused on evolution of 
spiral galaxies including gas (|Bottemall2QQ3l : iBaba et all 
[2009). There are also pure A/'-body simulations of 
disks with a lar ge number of particles such as 5 x 10^ 
(|Sellwoodll2010f ), but they still use two-dimensional ap- 
proximation and th e particle-mesh method w ith a grid 
size of 110 X 128 (jSellwood fc BinnevI 120021 ). Three- 
dimensional calculations focused on the evo lution of stel- 
lar bars (Athanassoula et al. 2005; Dubins ki et al.|[2QQ9l : 
ISellwood fc Debattistall2009f ). Thus, it is important to 
investigate the basic physics of pure stellar disks us- 
ing three-dimensional A/'-body simulations with a high 
enough resolution. 

In this paper, we report the result of high-resolution 
A^-body simulations of isolated stellar disks, in which 
multi-arm spirals spontaneously develop. We describe 
the method of our A/'-body simulations in Section 2. In 
Section 3, we show the results of simulations, and discuss 
the evolution of the spiral arms. We also discuss self- 
regulated mechanism of spiral arms and how the maxi- 
mum amplitude of spirals is determined. Section 4 is for 
summary and discussion. 

2. AT-BODY SIMULATIONS 

We performed a series of A^-body simulations of stellar 
disks in a fixed spherical dark halo potential. Our stellar 
disk a nd halo models are based on those in IBaba et"all 
(|2009[ ). Here, we briefiy summarize the parameters 
and how we generate the i nitial equilibrium dis ks. We 
adopted the NEW model (|Navarro et al.l I1997D as the 
dark halo model with the concentration parameter of the 
halo, c = 10. The virial radius of halo, i^h, is 122 kpc. 



and the mass within i?h, Mh is 6.4 x lO^^M©. ^ We 
adopted an exponential disk model as disk models. We 
varied total disk mass, M^, and initial Q at the refer- 
ence radius (8.6 kpc in our models), Qo- The scale ra- 
dius, i^d, is 3.4 kpc, and the scale height, Z(\^ is 0.34 kpc. 
We performed simulations with four different resolutions, 
A^ = 3xl0^,9xl0^ 3xl0^ lxlO^ and 3x10^ (hereafter, 
30M, 9M, 3M, IM, and 300k, respectively). We summa- 
rized our disk models in Table [H Hereafter, we regard 
model b (Md/Mh = 0.050 and Qo = 1-2) as our standard 
model. The circular velocity of the disk at = 10 kpc is 
about 200 km s~^ for model b. Its circular velocity pro- 
file is shown in Figure [TJ We g enerated initial d isk models 
using the Hernquist method (|Hernquistl 119931 ). Initially, 
the generated models are not exactly in an equilibrium 
at the galactic center. Thus, they cause ripples spreading 
through the disk from the center. To remove the ripples, 
we integrated the models for a few Gyrs randomizing az- 
imuthal positions of particle s every four steps to preven t 
the growth of spiral arms (|McMillan fc Dehnenl [20071) . 
After the ripples passed through the disk, we used it as 

the initial condition. 

We used a Burnes-Hut treeco de (iBarnes fc HuU [19861 : 
[Makino I 12004) on GRAPE -7 (iKawai et al.l I2006D and 
GRAPE-DR (Makino elini2007D. The opening angle, 
^, is 0.4 with the center-of-mass approximatio n. The 
maxi mum group size for a GRAPE calculation ([Makinol 
|1991[ ). ncrit, is 2048. For the time integration, we used 
a leapfrog integrator with a fixed stepsize of At = 0.29 
Myr for A' = 30M, 9M, and IM models and At = 0.15 
Myr for the other models. The gravitational potential 
is softened using the usual Plummer softening, with the 
softening length of e = 30 pc. This is small enough to 
resolve the typical spiral structures (~ 100 pc). For all 
runs, the energy error was less than 0.1% throughout 
the simulations. One may concerned about the drift of 
the disk in the halo potential because treecodes do not 
conserve the linear momentum of the system. Therefore, 
we investigated the position of the density center of the 
disk and confirmed that the drift does not occur (for the 
details, see Appendices A and B). 



TABLE 1 
Initial disk models 



Model 


Qo 


Md/Mh 


N 


a 


1.1 


0.050 


3M/300k 


b 


1.2 


0.050 


30M /9M /3M / IM /300k 


c 


1.3 


0.050 


3M/300k 


d 


1.4 


0.050 


3M/300k 


e 


1.5 


0.050 


3M/300k 


f 


1.8 


0.050 


3M/300k 


g 


0.5 


0.050 


3M 


h 


1.2 


0.075 


3M 


i 


1.2 


0.033 


3M 



We assume that a spherical region, where the mean density 
is 200 times as high as a background density, is virialized. As 
the background cosmology, we adopted a concordance cold dark 
matter (ACDM) model with parameters: Qm = 0-3, f^A = 0.7, 
and Hq = 70km s~^ Mpc"-*^. The formation redshift of the halo is 
set to be 1.0. 
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Fig. 1. — Circular velocity of model b (our standard model) as a 
function of the galactocentric distance. 

3. EVOLUTION OF STELLAR DISKS 

3.1. Evolution of spiral arms 

First, we show the evolution of spiral arms of our stan- 
dard model, model b. Figures [2] - [6] show the evolution 
of model b for AT = 30M, 9M, 3M, IM, and 300k, re- 
spectively. Top panels show the surface density of the 
disk in the Cartesian coordinate. Middle panels show 
the density contrast, in the polar coor- 

dinate. Here, T^{R^(j)) and t^{R) are averaged surface 
density in a polar grid from R to R -\- AR and from cf) 
to ^ + and a ring from R to R -\- AR, respectively. 
We used AR = 1 kpc and Acj) = 7r/64. Bottom panels 
show Fourier amplitudes at each radius obtained from a 
Fourier series: 

f^A^expHm^], (1) 

where m is the azimuthal wavenumber (i.e., the number 
of spiral arms) and Am is the Fourier amplitude. Only 
the amplitudes of m =2-6 are shown in this figure. Other 
modes are much smaller than these modes throughout 
the simulations. 

Initially many spiral arms with small amplitudes ap- 
pear (see t = 0.50 Gyr in Figures [2] -[6]). Eventually they 
merge and the amplitude of small wavenumber such as 
m = 4 becomes larger. Figure [7] shows the total power, 
the sum of squared amplitudes defined as X^m=i l^^l ^ 
at 7.5 ± 0.5 kpc. The data points are averaged over 
0.5 Gyr. In the beginning of the simulations, the total 
powers grow exponentially from their initial amplitudes, 
which are determined by Poisson noise. The growth 
timescale is ~ 0.4Gyr. Thi s value is compa rable to that 
in the case of the bar mode (|Dubinski et al.f ^OQ). In the 
case of N = 30M, we performed the simulation only to 5 
Gyr. However, the evolution in the run is quite similar 
to those in the runs with N = 3M, 9M, and 30M, al- 
though the total power grows from a smaller value. The 
total powers reach their peak values at t ~ 2 Gyr for 
N = 300k and t - 3 Gyr for N = IM. In the case of 
N — 3M, 9M, and 30M, the peak is not clear, but the to- 
tal powers have developed well at t > 6 Gyr. Clearly, the 
dependence on the number of particles exists. In the case 
of a larger number of particles, it takes more time for the 
spiral arms to develop because they start from a smaller 



amplitude of Poisson noise. After the spiral arms have 
developed, their number is consistent with that expe cted 
from the swing amplification theory (|Toomrelll981l ). in 
which spiral arms with 1 < kcvR/m < 2, where develop 
most efficiently. In our model b, m 4 for kcrR/m = 1.5. 

After t ^ 3 Gyr, the amplitudes of = 300k model 
start to decay (see Figure [7]), and the spiral arms almost 
disappear at t = 6 Gyr (see Figure [6]). The behavior of 
IM model is similar to the 300k model. On the other 
hand, the spiral arms in N =3M, 9M, and 30M models 
are still prominent after 6 Gyr, and their amplitudes are 
\Am\ ~ 0.05 even at the end of the simulation, t = 10 
Gyr, in contrast to \Am\ ^ 0.03 in A" = 300k model (see 
Figure |4]) 

Evolutions of the amplitudes of spiral arms are quali- 
tatively different between higher and low resolutions. In 
A^ = 300k model, the amplitude grows rapidly in the 
first 2 Gyr, and after this rapid growth it decreases fairly 
rapidly. On the other hand, the amplitude in A" = 3M, 
9M, and 30M models keeps growing to the end of simu- 
lation, i.e., t = 10 Gyr. The highest resolution run may 
show the same evolution in > 5 Gyr. 

We found that spiral arms developed in the disk are not 
stationary, but timedependent. As shown in Figure [71 
the amplitude of the arms in A" = 3M and 9M models at 
R = 7.5 kpc oscillates quasi-periodically in the timescale 
of ~ 1 Gyr. Figures [2] - [6] show that all modes of spiral 
arms are timedependent and the dominant mode changes 
spatially. Each spiral arm is wound up because of the 
galactic differential rotation. As a result, the global spi- 
ral arms break up into smaller fragments whose sizes are 
typically a few kpcs. These fragments eventually collide 
and reconnect with other fragments due to the differen- 
tial rotation, and the global spiral arms revive. This pro- 
cess of breaking up and reconnection repeats throughout 
the simulations. 

Figure [8] shows the radial distribution of the surface 
density, S, radial velocity dispersion, o-r, Toomre's Q, 
and scale height, (2:^)^/^, of model b for A" = 3M and 
300k. These values are averaged over bins of 500 pc. 
While the distribution of and Q at t = 10 Gyr are very 
similar in both models, (z^)^/^ in A" = 300k increases 
more rapidly than that in A^ = 3M. The evolution of 
(2:^)^/^ is caused by the two-body relaxation. The two- 
body relaxation time of this model is ~ 3 Gyr and ~ 30 
Gyr for A" = 300k and 3M. The evolutions of aR and Q 
are also faster in A^ = 300k models. We will discuss the 
effects of the number of particles in Section 3.5. 

3.2. Evolution of the Q value due to the spiral heating 

In this section, we compare the results of models with 
different initial values of Q, Qo, and investigate how the 
Q evolve in time. We performed simulations for models 
with different Qo (model a-f). Top panels of Figure [9] 
show the time (a) evolution of Q, (b) the growth rate 
of Q, i.e., dQ/dt, and (c) the total power of the modes, 
which corresponds to the amplitude of spiral 
arms. These values of Q and dQ/dt are averaged over 
the range of 5-10 kpc and over 0.5 Gyr. The total powers 
are evaluated at 7.5 ±0.5 kpc. If the initial disk is colder 
(i.e., smaller Qo), the amplitudes of spiral arms tend to 
be larger in both A" = 300k and 3M models. It is also 
clearly seen that Q increases more rapidly in colder initial 
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Fig. 2. — Evolution of spiral arms for model b with N = 30M. Top panels show the surface density, 
density normalized at each radius, and bottom panels show the Fourier amplitudes. 
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disks (i.e., smaller Qo). In all N = 3M models, the 
amplitudes of spirals tend to increase toward t = 10 Gyr, 
except for model a. On the other hand, the amplitudes 
start to decrease soon after the simulations start in all 
N = 300k models (right bottom panel of Figure [9]). The 
peak amplitude is larger for models with smaller Qo, for 
both numbers of particles. 

From comparison between panels (a) and (c) of Figure 
[9l it seems that the increase of the Q is caused by the 
spiral arms, since the Q rapidly increases when the am- 
plitude of spiral arms is large. Since the surface density, 
S, and epicycle frequency, a^, do not change significantly 
throughout the simulations (see Figures [8]), the change 
of Q depends only on the radial velocity dispersion, aR, 
by the definition of Q: 



3.36GS' 



(2) 



Therefore, dQ/dt can be interpreted as a "heating rate." 
The evolution of dQ/dt and that of the total power of the 
spiral arms, ^ l^mP, are very similar (see the left panels 
(b) and (c) of Figure [9]), indicating that the spiral arms 
increase the velocity dispersion of stars. The similarity 
between dQ/dt and ^ |^mP is also visible in the 300k 
models (right panels (b) and (c) of Figure [9j) . Therefore, 
this mechanism seems to be independent of the number 



of particles. 

In order to confirm the hypothesis that the amplitude 
of spiral arms, determines the heating rate, 

dQ/dt, we analytically estimate dQ/dt from the ampli- 
tude of spiral arms in the simulations. The relation that 
dQ / dt is proportional to the square o f amplitudes is sug- 
gested bv lCarlberg fc Sellwoodl (|1985f ) . They derived this 
relation by considering the perturbing potential of spiral 
arms. We derive the relation between dQ/dt and the 
spiral amplitude in a different way. As shown in the 
previous section, the global spiral arms are transient; 
splitting to smaller sub-arms and merging into global 
arms recurrently occur. Therefore, as a first approxi- 
mation, they behave like 'material arms' consisting of 
several massive clumps. With this assumption, we can 
estimate the increase of stellar velocity dispersion using 
the same equation as that describes dy namical heatin g 
of stars by GMCs in galactic disks (K okubo fc Ida|[T992l ), 
replacing a spiral arm in our simulation with several mas- 
sive clumps. As discussed above, the evolution of Q cor- 
responds to that of the velocity dispersion. We there- 
fore discuss the time derivative of the velocity dispersion 
below. The relaxation time of disk stars due to the dy- 
namical heating by clumps with mass Mc, whose number 
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Fig. 3.— Same as Figure E] but for Af = 9M. 



density is ric, is given by 



7rncG2(Mc + ms)2lnA' 



(3) 



where v and mg are the three-dimensional velocity disper- 
sion of disk stars and mass of a star respectively, a nd In A 
is Coulomb logarithm (jBinnev fc Tremaiii3l2QQ8f ). Here, 
since Mc rris, we can neglect mg. The number density 
of stars in spiral arms is given by ric = Sc/((^^)"'^/^Mc), 
where (z^)^/^ is the scale height of the disk and Sc is 
the surface density of the spiral arms. The scale height 
^^2\^i/2 -g giygj^ l^z^Y^^ ^ CFz/^i where cFz is the ver- 
tical velocity dispersion and Vt is the a ngular speed of 
the d isk. In a disk system, cjz — v (e.g., iKokubo fc Idal 
Il992f ). the relaxation time of a disk can be written as 
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clump of the spiral arms are given by 



m 

= S ^ ^ \Aj^ 



(6) 
(7) 



where and S are the total mass and surface density of 
the disk, and m and Ara are the number and amplitude 
of the spiral arms. Substituting Equations (jl]), (jS]), and 
([7j) to Equation ([5j), we obtain 



dv 
'dt 



TrSG^Md^^lnAg \A. 



2v^ 



TrScG^Mcl^lnA' 

By definition, the relaxation time is tg = ^^^y^- 
change of velocity dispersion is written as 

dv 1 dv'^ V 
'dt ~ 2v~dt ~ 2U' 



(4) 
The 

(5) 



From Equation (j2j), the time derivative of Q is 

dQ K, d(jR 

~dt 



(8) 



(9) 



3.36GS dt ' 

where we assumed that n and S are constants. As- 
suming that the three-dimensional velocity dispersion is 
V = \/?t(jR^ from Equations (|8]) and (|9]), we obtain 

(10) 



We assume that the mass and surface density of each 



dQ ^ 
dt 

If we scale as R 



8(kpc) = 1, Md = 



3.2 X lO'^M© = 1, 
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and Vt — 190km s 

^^3.91nA 

dt 



^ = 1, then 
190 [km s-^] 



^^Gyr-i(ll) 



where G = 0.39 and c:^ 1.5 in our model. Using Equa- 
tion (pT]) and the amphtudes, obtained from the 
simulation, we can estimate dQ/dt. We take the sum of 
Fourier components with m =4-6, which are dominant 
modes. As can be seen in the snapshots and maps of 
T>{(j)^ R)/T>{R) (Figures [2]-[6]), it would be unphysical to 
include modes with m < 3. We adopted In A = 1.0 be- 
cause the scale height of the disk is comparable to the 
size of the clumps. Figures [10] and [TT] show the compar- 
ison between the analytic and the numerical results of 
model a and b. It is clear that behavior of dQ/dt in the 
simulations is quantitatively reproduced by the analytic 
estimate in models a and b, for both N = 3M and 300k. 
Thus, we conclude that scattering of stars by spiral arms 
can heat up the stellar disks, and that the heating rate is 
proportional to the squared amplitude of the spiral arms. 

3.3. The maximum amplitude of spiral arms 

In the previous section, we showed that dQ/dt is tightly 
coupled with the amplitude of spiral arms, and this cou- 
pling is well understood in a simple picture that spiral 



arms gravitationally scatter stars and increase Q. What 
controls the amplitude of spiral arms? Figure [9] shows 
that the amplitudes of spiral arms decay, as Q increases. 
This suggests that the amplitude of spiral arms is deter- 
mined by Q. 

Figure [12] shows the evolution of models in the plane 
of Q and the total power of the spiral arms. These val- 
ues are measured at 7.5 kpc and averaged over 0.5 Gyr 
in the same way as the previous section. At the be- 
ginning of the simulation, both Q and the total power 
are small. The amplitude of spiral arms grows rapidly, 
and the models move in the right-upward direction as 
shown by arrow 1. Once the amplitudes reach their peak 
values, they decay and models move right- downward (ar- 
row 2). In this phase, the trajectory in the Q-Y^ |^mP 
plane seems to follow a roughly straight line, irrespective 
of models and the number of particles. In other words, 
there seems to be a "forbidden region" in the left-top of 
the Q-Y^ l^mP plane, where both amplitudes and Q are 
large. This result implies that the maximum amplitude 
is determined by Q. At the beginning, the amplitude is 
smaller than the maximum amplitude and therefore the 
spiral arms can grow with time. However, once the am- 
plitude reaches to its limit, it starts to decay because Q 
increases due to heating by spiral arms, and the maxi- 
mum amplitude decreases. 
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Here we try to estimate the maximum amplitude of 
spiral arms, assuming that the spiral arms grow through 
the collapse due to the gravitational instability, and that 
they evolve until they reach an approximate dynamical 
equilibrium. Under this assumption, the amplitudes can 
be simply estimated as a density contrast before and after 
the collapse. We assume that stars in a region of the 
disk with a size of a critical wavelength, Acr, collapse to 
a spiral arm. 

The initial energy of the region to collapse can be ex- 
pressed as 



Eo = Ko^Wo = -Ma^ 



C 



(12) 



where M, a, and ro are the mass, velocity dispersion, and 
radius of the system, and C is a fixed value. We treat 
C as a parameter, depending on geometry and density 
distribution; e.g., for a homogeneous sphere, C = 3/5. 
Assuming that the collapsed region is virialized when it 
forms a spiral arm, the arm satisfies the virial theorem 
and the potential energy after the virialization is 



(13) 



The amplitude of the spiral arms is obtained from the 
density contrast of the initial and virialized density. If 
we assume that the virialized density is the mean density 



inside the half-mass radius, rh, we obtain the amplitude 
from the ratio of the initial and virialized densities, 



p _ M/(2rg) 



1 



M/rl 2 (0.45GM2)3 



|2(i^o + Wo)|', (14) 



where we adopt the half-mass radius, rh = 0.45GM^/| VK| 
(Binney & Tremaine 2008), and W is obtained from 
Equation (p!3|) . From Equation ([2]) and the critical wave- 
length, Acr = 47r^Gi;//T:^, we obtain 



4 



(3.36)^ 
47r2 



GSAcrQ • 



(15) 



Assuming that a = \/2>(Jr and the radius of the sphere, 
''o = Acr/2, we can rewrite Equation (fTl)) using Equations 
(jfe)) and ^ as 



^ = i('4.4C-0.95^g2 



(16) 



where we assumed that Eq < 0. Since M ~ SA^^, the 
density contrast is written as a function of Q, 



= - (4.4C-0.95g2)^ 



(17) 



po 2 

While the density contrast in Equation (fTTjl is defined in 
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Fig. 7. — Left panel: time evolution of total power 

(Em=i l^ml^) for model b at = 7.5 kpc. Right panel: same as the left panel, but 
in the logarithmic scale. The solid thin line shows an exponential timescale of 0.37 Gyr. 
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Fig. 8. — Evolution of surface density, E, radial velocity dispersion, ctr^ Toomre's Q value, and scale height, (z^)-*^/^, of disk for Model 
al, = 3M (left) and 300k (right). Here, we defined the rms of the vertical position of stars as l^z^)^l^ . 
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Fig. 9. — Time evolution of Q (a), dQ/dt (b), and total power (c) for model b, with = 3M (left) and 300k (right). These Q and dQ/dt 
are averaged between 5-10 kpc, and the total powers are evaluated at 7.5 ± 0.5 kpc. 
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Fig. 10. — Comparison of dQ/dt between that obtained from simulations (measured) and that estimated from Equation (|11|) (analytic) 
for model a, = 3M (left) and = 300k (right) at = 7.5 kpc. 
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three dimensions, the amphtude we obtained from the 
simulation is one-dimensional because it is the contrast 
of the radially averaged surface density. Therefore, we 
define the amplitude of a spiral arm as 



An 



1/3 



1. 



(18) 



From Equations (pT|) and (fT8|) , the amplitude relates with 
Q as 



3.5C- 1.0-0.75Q^ 



(19) 



The black curve in Figure [T2l shows (O.IA^)^ obtained 
from Equation ([19]), where we assumed C = 1.0. Al- 
though Equation ([T9]) qualitatively explains the ampli- 
tude as a function of Q after it reaches the maximum, 
the amplitude obtained from the simulations is smaller 
than that obtained from Equation (fT9|) by a factor of 10. 
The possible reasons are as follows. (1) We assumed that 
a homogeneous region collapses to estimate the density 
contrast, but this is not the case in a disk. Especially, the 
scale height of the disk is much smaller than ro , which we 
assumed as the initial radius. (2) We used the radially av- 
eraged surface density to calculate the Fourier amplitude 
from the simulations. This treatment may underestimate 
the local amplitude of spiral arms because the averaged 
density depends on the radial width for averaging. We 
chose a radial width smaller than the critical wavelength. 
However, the amplitudes increased by ~ 10% when we 
halved the width. (3) Growth of spiral arms does not 
complete because the galactic shear breaks up the spiral 
arms before they are virialized. 




Q 

Fig. 12. — Relation between Q and total power of spiral ampli- 
tude. The value is evaluated at 7.5 kpc and averaged in 0.5 Gyr. 
Black curve shows (O.lAm)^ obtained from equation (|19|) . Arrows 
shows the direction of the time evolution. 

We do not insist that our simple model gives a fully 
correct description of the mechanism through which Q 
controls the amplitude of spiral arms, but it is clear from 
Figure [12] that Q determines the amplitude. Thus, the 
spiral arms evolve in a self-regulated manner as follows. 
Spiral arms grow fron i small density perturbations by the 
swing amplification (jToomrd Il981f ) to their maximum 
amplitudes limited by Q. The spiral arms scatter disk 
stars, and as a result the velocity dispersion of the disk 
star increases. The heating rate, dQ/dt, is proportional 



to the squared amplitudes of spiral arms (see Equation 
(fTT]) ). This heating mechanism increases Q, and there- 
fore the amplitude of spiral arms decreases. As a result, 
the heating rate decreases as Q increases. Through this 
evolution, the spiral arms become asymptotically faint as 
qualitatively shown by the black line in Figure [121 but 
its timescale is comparable to the cosmological time, i.e., 
10 Gyr. 

3.4. Effects of the initial Q value and disk mass 

In order to see how the initial Q, and disk mass fraction 
affect the evolution and morphology of spiral arms, we 
performed three additional runs: an unstable disk (Qo = 
0.5; model g), a massive disk (Md/Mh = 0.075, model 
h), and a less massive disk (M^/Mh = 0.03, model i). In 
all models, the number of particles is 3 x 10^. Figure [13] 
shows the snapshots of models f-i. Model f is a model 
with a large initial Q, Qo = 1-8 and shown in Section 
3.2. We show it again as an example of a model with a 
large Qq. 

Model g is initially cold and unstable, therefore strong 
spiral arms develop in the first 0.5 Gyr. As shown in 
top panel of Figure [HI the disk is soon heated up to 
Q ~ 1.6, and the Q keeps increasing. As expected from 
our theory in Section 3.3, the amplitude of spirals then 
decreases quickly, and they are very weak at t = 6 Gyr. 
The final density and velocity profiles are quite different 
from original ones. In model f, whose disk is initially hot 
{Qo = 1.8), spiral arms do not develop, and therefore Q 
(or the velocity dispersion) stays nearly constant (see the 
left panels of Figure [9]) . 

Models h and i have the same parameters as those 
of model b (standard mo del) except the disk mass ra - 
tio to halo, Md/Mh . iCarlberg fc Freedmanl (|1985D 
showed that the number of spiral arms in numerical sim- 
ulation is consistent with that pre dicted by the swing 
amplification theory (Too mrd Il981f ). and massive disks 
have a smaller number of spiral arms. Our results 
are consistent with this previous result. We can es- 
timate the number of spiral arms as follows. The 
swing amplification is characterized by a parameter, 
X = kcrR/m = t<?R/2T:GT^m, where k^r is the critical 
wavenumber (jBinney fc Tremaine 2008). Spiral arms de - 
velop most effectively when 1 < X < 2 (|Toomr jfigSlf ). 
Therefore, we can estimate the dominating number of 
spiral arms, m, as 



^^R 



27rGSX 47rGE ' 



(20) 



where we adopted X 2. In the case of model b, we 
obtain m = 6 at 8 kpc, which roughly agrees with the 
result of the simulation. This estimate is applicable for 
other models with different disk mass fractions. Since the 
halo mass, Mh, is fixed in our models, the surface density, 
S, is proportional to the disk mass fraction, Md/Mh. The 
numbers of spiral arms of models h and i are estimated 
as m = 4 and m = 9 at 8 kpc, and they also agree with 
the results of the simulations. 

Figure [14] shows the time evolution of Q averaged 
within 5-10 kpc (top) and the total power at = 7.5 
kpc (bottom), respectively. As was the case with mod- 
els a-f, the large total powers correspond to the rapid 
increase of Q. We investigated the time evolutions of 
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dQ/dt from the simulations and evaluated them using 
Equation ([9]) in the same way as models a-f. Figure [T5l 
shows the results. The analytic results again agree well 
with the simulations in these models. 

In fact, the evolution of these models are qualitatively 
similar, but different in details. The evolution of model 
h, which has a more massive disk, is faster than that 
of model b (our standard model). In model h, the total 
power of the spiral arms grows to ^ |^mP ^ 0.02, at t ~ 
2-4 Gyr, whereas it is 3 x 10~^ in model b. This difference 
causes the faster decay of the amplitude in model h (see 
bottom panel of Figure [T4|). Furthermore, the number of 
spiral arms decreases in model h. Initially, the number 
is around four (see 2.0 Gyr in Figure [T3|). but three at 
6.0 Gyr. In model h, the effective angular-momentum 
transport occurs due to asymmetric structures in the disk 
and the surface density of the inner region increases. This 
reduces the number of the spiral arms. 

3.5. Effects of the number of particles 



As we showed in Figures [2] - [6l while the spiral arms 
survive for more than 10 Gyr in = 3M, 9M, and 30M 
models, this is not the case in = 300k and IM models. 
Although the disks in A' = 300k models become feature- 
less after t = 6 Gyr (Figure [6]), this is not mainly by 
the effect of the two-body relaxation. We can estimate 
the heating rate due to the two-body relaxation from 
the result of model f, where spiral arms do not develop, 
the relative heating rate for model f with A^ = 300k at 
Q = 1.8 is around 0.5%/Gyr. For a disk, the relative 
heating rate is proportional to Q~^. Thus, for Q = 1.2 
and A" = 300k, the heating rate due to the two-body 
relaxation is 3%/Gyr, which is much smaller than the 
actual increase obtained from simulations. Heating due 
to spiral arms is dominant even in A" = 300k models. 

The major differences between A" =3M-30M models 
and A^ =300k and IM models are (1) the time when the 
spiral arms reach their peak amplitudes and (2) the peak 
amplitude themselves. The spiral arms initially grow 
from small perturbations of the density originated in the 
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Fig. 14. — Time evolution of Q averaged within 5-10 kpc (top) 
and total power, X]m=i l^rn|^, averaged in 1 kpc at 7.5 kpc (bot- 
tom) for models b, g, h, and i, at = 7.5 (kpc). 

Poisson noise of disk stars and grow up to their maxi- 
mum amplitude determined by Q. For the same initial 

N = 300k models reach the maximum amplitudes 
much faster, and the values are large. A smaller num- 
ber of particles generates a larger Poisson noise. There- 
fore, spiral arms in N = 300k and IM models can reach 
its maximum amplitude fast er than those in N = 3M- 
30M models (jSellwoodll20To [). This means that the peak 
amplitude in N = 300k models is larger than those in 

= 3M-30M models, because they can develop before 
the Q becomes large. The larger peak amplitude results 
in larger heating rate. Thus, the Q of the N = 300k 
models increases faster than that of the N = 3M-30M 
models. 

As shown in Figure [71 the evolution of the models with 
N =9M and 3M is not significantly different after t ^ 2 
Gyr. In other words, once the amplitude of the spirals 
reaches a certain level P - 2 X 10-^), the differ- 

ences in the initial conditions are no longer important. 
In section 3.1, we estimate the timescale of the initial 
exponential growth as around 0.4 Gyr. If we start from 
the noise level of real galaxies, ^ 10^^ particles, it will 
take 1.6 Gyr longer than in the case of A^ = 9M. This is 
still shorter than the cosmological timescale. 

The two-body relaxation has serious effects on the 
heating of the disk only when the number of particle is 
quite small. For example, the number of particles used in 
the simulation of .Sellwood fc Carlberg (1984) was only 



2 X 10^. Since the heating timescale of our A" = 300k 
model is ~ 10 Gyr for Q = 1.0, the relaxation time of 
their model would be ^ 1 Gyr. Thus, it seems that in 
ISellwood &: Carlberd ([1981 spiral arms are weakened by 
the heating due to the two-body relaxation. 

4. SUMMARY AND DISCUSSION 

4.1. Summary 

We performed three-dimensional A'-body simulations 
of pure stellar disks with spiral arms and investigated 
their dynamical evolution. We confirmed that the spiral 
arms are transient and recurrent. Contrary to previous 
results, we found that spiral arms in pure stellar disks 
can survive for more than 10 Gyrs, when we use a suf- 
ficiently large number of particles. We also found that 
spiral arms of a stellar disk are self-regulated. The spiral 
arms grow by the swing amplification up to their maxi- 
mum amplitudes determined by Toomre's Q value at the 
moment. The amplitude becomes smaller as Q increases. 
The spiral arms heat up the disk, or increase the veloc- 
ity dispersion of stars, by scattering the disk stars. As a 
result, Q increases, and the amplitude of spirals is sup- 
pressed. We found that the heating rate, which is given 
by dQ/dt^ is roughly proportional to the square of am- 
plitudes. It means that the heating timescale becomes 
longer as Q increases. Thus, the spiral arms heat up the 
stellar disk and increases Q, but at the same time the in- 
crease of Q results in the decay of the spiral amplitudes 
and a smaller heating rate. This self-regulating relation 
among Q, spiral amplitudes, and the heating rate main- 
tains the spiral arms for more than 10 Gyr. 

In the case of the smaller number of particles (A" = 
300k), however, the spiral arms become faint much faster 
than in the model with A^ = 3M, 9M, and 30M. We found 
that the initial exponential growth of the density pertur- 
bation depends on the number of particles (see Figure 
[7j) . Spiral arms initially grow from density perturbations 
originated from the Poisson noise of the initial condition 
through the swing amplification. A smaller number of 
particles results in a larger seed noise. Therefore, spiral 
arms grow faster up to their maximum amplitude deter- 
mined by Q. The rapid growth of the amplitude causes 
a rapid heating of the disk and also rapid decay of the 
spiral arms. 

Our results show that the timescale of the initial ex- 
ponential growth is 0.4 Gyr. It means that even if we 
start with smooth disks A^ ~ 10^^ stars, the total power 
of arms reach the level of ~ 2 x 10~^ in a few Gyr. In real 
galaxies, moreover, the disks may be initially perturbed 
by hierarchical mergers and/or GMCs. If the initial per- 
turbation is comparable to or larger than this level, the 
spiral arms that look very similar to those in real spiral 
galaxies are formed within a Gyr or less. It is beyond 
the scope of the present paper to discuss the amplitude 
and shape of the initial perturbation in a realistic situa- 
tion, which should be ultimately studied in a cosmologi- 
cal context with a sufficiently high numerical resolution, 
e.g., the generated disk should be represented by larger 
number of particles to avoid numerical artifacts. 

4.2. Effects of Gas and Star Formation 

In this paper, we showed that the presence of gas is not 
essential in maintaining spiral arms. However, real spi- 
ral galaxies have gas and its effect is not negligible. Gas 
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5. — Comparison of dQ/dt between measured and analytic results for models g (left) and h (right), = 3M R = 7.5 ± 0.5 kpc. 



Fig. 1 

in galactic disks can work as both cooling and heating 
sources. The velocity dispersion of gas is smaller than 
that of stars because gas clumps lose their kinetic energy 
due to dissipation when t hey co llide each other. As sug- 
gested bv .Sellwood fc CarlbergI ((1984), new stars formed 
from the gas have smaller velocity dispersions than that 
of old disk stars because of the smaller velocity disper- 
sion of gas. The new stars can dynamically cool the 
disk. Moreover, the existence of gas reduces the effec- 
tive Q of the disk (Jog & Solomon 1984). On the other 
hand, the smaller Q ma kes the disk dynamically unstable 
fBertin fc Romeolll988[ ) and it may cause faster heating 
of the disk. In addition, gas trapped into stellar spiral 
arms would cool due to the dissipation and its gravity 
may strengthen the amplitude of spiral arms. The larger 
amplitude of spiral arms will cause faster heating of the 
disk. However, even if the gas enhances the amplitude of 



spiral arms, the self-regulating mechanism that we sug- 
gested in this paper will determine the spiral amplitudes. 
Therefore, the disk would keep the spiral arms for a long 
time as in the case of pure stellar disks. In reality, the 
interaction between the gas and stellar spirals is com- 
plicated. We are now performing simulations with gas 
and will present the results elsewhere (Baba et al. in 
preparation). 
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APPENDIX 

EFFECTS OF OPENING ANGLE OF TREECODE 

We used a treecode for calculation of forces. Since a treecode does not conserve the linear momentum of systems, the 
disks might drift during the simulations. Since we used a fixed background potential as a halo model, the drift from 
the center of the halo might result in an artificial m = 1 perturbation. This effect depends on the opening angle of 
treecodes, ^, which is a parameter in the tree approximation. Smaller can achieve more computational accuracy, but 
requires longer CPU time. If the drift is significant to induce the artificial perturbation, we expect larger amplitude of 
spirals developing in the later phase for larger 6. In order to confirm this, we performed simulations with 6 = 0.2, 0.4, 
and 0.75 for model b with N = 300k and investigated their convergence. Figure [16] shows the evolution of the total 
powers. They do not show any clear differences. We also confirmed that the amplitudes of the m = 1 mode do not 
grow throughout the simulations (see Figure [T7|) . In the results in Section 3, we adopted = 0.4 for all simulations. 

THE TIME EVOLUTION OF THE FOURIER AMPLITUDES AT R < 1 KPC 

In this section, we show the time evolution of the Fourier amplitudes at < 1 kpc of the disk, which corresponds 
to the deviation of the density center of the disk from the center of the halo potential, in order to confirm that the 
drift of the density center of the disk is not significant. 

We investigated the deviation of the density center of the disk from the origin, which is the center of the halo 
potential, on the x-y plane. This directly shows the drift of the disk. The dotted curve in Figure [18] shows the distance 
of the density center from the origin at R < 1 kpc for model b with = 30 0k, Ar. The density center is calculated 
from the highest local density using the method of iCasertano fc HutI (|1985 t). Figure [TS] shows that the deviation of 
the density center is as small as the softening length of 0.03 kpc and does not increase during the simulation. The 
full curve in Figure [18] shows the time evolution of the Fourier amplitudes of m = 1, Ai. It shows that Ai traces the 
deviation of the density center, Ar. The amplitude is due to the Poisson noise of the particle distribution. Hereafter, 
therefore, we adopt Ai at < 1 kpc as the index of the deviation of the disk center from the halo center. 

Figure [19] shows the Fourier amplitudes of m =1-6 at < 1 kpc for model b. These amplitudes are averaged over 
0.5 Gyr. For comparison, we also plotted the result of model f with N = 3M, in which spiral arms do not develop, 
because the initial Q is large. Initially, all amplitudes show the fiuct nations due to Poisson noise. In the case of small 
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t (Gyr) 

Fig. 16. — Evolutions of total powers dX R = 7.5 ±0.5 kpc for model b with = 300k. Solid, dashed, and dotted curves show the results 
with opening angles, 6 = 0.2,0.4, and 0.75, respectively. 




t (Gyr) 

Fig. 17. — Time evolution of the amplitude of m = 1, Ai, at = 7.5 ± 0.5 kpc for model b with = 300k. Solid, dashed, and dotted 
curves show the results with opening angles, 9 = 0.2,0.4, and 0.75, respectively. 

numbers of particles {N < IM), the amplitudes do not increase during the simulations. In the case of larger numbers 
of particles {N > 3M), on the other hand, the amplitudes of m = 1-3 modes increase after the spiral arms developed 
(after ~ 2 Gyr). We consider that this increase of the amplitudes at < 1 kpc is caused by the spiral arms developed 
in the outer part of the disk (e.g., at 7.5 kpc), because the amplitudes does not increase in the case of the disk without 
spiral arms (see model f in figure [19]). Moreover, it is unlikely that these small increases of the amplitude at R < 1 
kpc affect the evolution of the spiral arms in the outer part of the disk. In model b with N = 30M, for example, the 
amplitude grows to ^ 30 times as much as the initial fiuctuation during the first 2 Gyr (see figure [7j). After 2 Gyr, 
the amplitudes at < 1 kpc start to increase. If the increase of the amplitudes at < 1 kpc cause the evolution of 
the spiral arms of the outer part, the increase during 2-4 Gyr at R = 7.5 kpc must be larger than those of the first 
2 Gyr. In fact, however, the spiral amplitudes at = 7.5 kpc grow to only 2-3 times during 2-4 Gyr. Therefore, we 
can safely conclude that the increase of the amplitudes at < 1 kpc is caused by the asymmetry of the self-developed 
spiral arms in the disk. 
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Fig. 18. — Time evolution of the amplitude of m ■ 
model b with = 300k and 6 = 0.4. 



1, Ai, and the distance of the density center from the origin at < 1 kpc, Ar, for 
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Fig. 19. — Time evolution of the Fourier amplitudes at < 1 kpc. 



